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First results on the autocorrelation behaviour of a recently proposed fermion algorithm by M. Liischer are 
presented and discussed. The occurence of unexpected large autocorrelation times is explained. Possible improve- 
ments are discussed. 
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1. INTRODUCTION 

The basic ideas of the method are described 
in [1] and in this volume in [2] and will not be 
discussed here. This paper uses the notation of 
[2]. The tests were done in a 423 SU(2) local 
gauge theory with 2 flavours of Wilson quarks 
with periodic boundary conditions. 

2. CURRENT IMPLEMENTATION 

Recall that the action of the approximating 
theory is given by 



Sk[U,4>] = Sg[U] + 

with 

Q = fsiD + m)/[cM{8 + m)], 



(1) 



(2) 



and D the Wilson-Dirac operator, cm > 1 is an 
arbitrary constant to avoid large bosonic autocor- 
relation times. 

The local structure of the interactions allows 
us to use conventional MC procedures for gauge 
and bosonic fields. Namely writing (neglecting 
constant terms) 

Sb[4>k,x\ = 4'k,x Ak4>k,x + Bk,x 4>k,x + h.C. (3) 

with Ak constant and diagonal (in the chiral rep- 
resentation), one can use the standard gaussian 
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over-relaxation and heatbath methods. For the 
gauge fields one may write 



S[U,,.] = ti \U,,. ( -^F, + F,,,) \, (4) 



where Fg are the staples of the gauge action and 
i^^^a; is the induced action of the bosonic fields 
(which is again a quaternion, i.e. can be written 
as a coefficient times an SU(2) matrix). Writ- 
ten this way, one can again apply heatbath [3,4] 
and the standard overrelaxation methods for the 
SU{2) (and generally for SU{N)) case. 

In our implementation, an iteration is made up 
of one bosonic heatbath sweep, N^^ bosonic over- 
relaxation sweeps, one gauge heatbath sweep and 
finally N^^ gauge overrelaxation sweeps. Recall 
the gauge case, where TV^j. can be used to decrease 
Tint and the dynamical critical exponents [5] sub- 
stantially. 

The implementation was done on a Quadrics 
Q16 machine with 128 nodes. On one node, the 
algorithm needs for the bosonic over-relaxation 
step 170/xs per update, field and site and 250/xs 
for the heatbath and gauge over-relaxation steps. 
This allowed us to accumulate up to 200K iter- 
ations for each data point. The maximum lat- 
tice size which can be reasonably simulated (us- 
ing preconditioning) on such a machine is roughly 
a 16'* lattice with 100 bosonic fields. On this lat- 
tice a sweep would need approximately 5s with a 
memory consumption of roughly 115 MW. 
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Using dQ / d\\\=i = Q — C075, we write 
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Figure 1. rmt against the number of fields for a 
4^ X 8 lattice. The bare parameters are /3 = 1.75 
and K = 0.165. N^,. is 1 in all cases. 



3. AUTOCORRELATION TIMES 

We performed measurements on 4^ x 8 and 
6^ X 12 lattices, holding roughly TO,r/"^p ~ 0.93 
and rriTrL « 7. Figure 1 shows the integrated 
autocorrelation time against the number of the 
bosonic fields n, holding all other parameters 
fixed. One sees clearly two points: first, rmt 
is roughly proportional to n and secondly, over- 
relaxing the gauge field does not decrease rmt (it 
even increases rmt in units of CPU time linearly). 
Furthermore, the proportionality factor is rather 
high. 

3.1. n-behaviour of Tint 

Recall that 



< A,B >= -ti {AB'^} 



(5) 



defines a scalar product on the space of quater- 
nions, so that we may talk about angles and vec- 
tors. Let us look at the expectation value of 
tr{t/'^,a;-F^_j,}. We write here Q as co75{l - kH} 
and replace t/^^a; -^ ^Uij,,x with A G R and write 



E.,,tr{t/,,.i^^,J 



SA" 



lA = l 



32(co/c)2 ^^ 



(6) 



The last term comes from the X'^U^^xUJ^ x expres- 
sions in Q^ and, as we will see, causes the trouble. 
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SA" 



.) 



(7) 



tr {^^^(Q(Q - C075) + (Q - co75)Q)}) . 

P'{x)/ P{x) approximates —1/x for (n, e) — ^ 
(00, 0) for all polynomials we use, so we can write 



2V 



tr 



{T} 



B(eoK;)^ 



Z-^k 



^l<i>k 



(8) 



The first two terms are constant in n, but the 
third one is a sum of n nearly independent vari- 
ables and thus rises linearly with n (this can be 
explicitly calculated in the free case). This means 
that ||-F^,a;|| rises linearly with n, which can in- 
deed be observed in the simulation. This explains 
the unwanted behaviour: after a bosonic update 
the force induced by the (j)}. fields points roughly 
in the same direction as the gauge field, with a 
length proportional to n and an angle of order 
l/\/n. The following gauge update allows the 
gauge field to move only by an angle of the same 
order. Holding all parameters except n fixed, we 
may assume, that rmt is proportional to the time 
needed for the gauge field to turn an angle of or- 
der 1. Since we perform a random walk this is 
proportional to n, which is exactly what we see. 
Also multiple updating of the gauge link does not 
improve the situation since the gauge link will 
only fluctuate narrowly around the large force, 
therefore N^^ has indeed no effect. This situation 
seems to persist as long as we update the gauge 
and bosonic degrees of freedom independently. 

4. IMPROVEMENTS 

There are two immediate ideas on how to im- 
prove the situation. Firstly, one has to keep n as 
small as possible and secondly, updating gauge 
and bosonic fields together may improve the be- 
haviour. The first possibility is already discussed 
in [2], so we will present an attempt on the second 
point which is presently being tested. 

4.1. Modified Updating 

One obvious possibility to increase the freedom 
of the gauge links is to update the link and some 
bosonic fields together, so that U^^^x and i^^^j, turn 




Figure 2. Distribution of the angle between up- 
dated and original gauge links using heatbath up- 
dating for different n. The leftmost curves are 
done using standard update, the rightmost ones 
using a combined update. 



simultaneously. One way to achieve this is to in- 
tegrate out the bosonic fields at the endpoints of 
the gauge link and update the link according to 
the effective action 



e-5,„,,(!7^,.) ^ 

/ Ilk <i4>k,x<i4>k,x+fj,e~^'^^'' 







(9) 



via the heatbath procedure. Since the integrals 
in 4> are gaussian, one can show that 



SeS,g{U^,x) = tr {u^,xFl,} , 



(10) 



with i^^^a; of the same form as -F^,a;, but slightly 
more complicated. To fulfill the stability equa- 
tion Jdud4>p{u, 4>)T{u, 4> -^ u', 4>') = p{u', 4>') , 

one has to first do a heatbath update t/^^a; ~^5,ff 
U' J,, then update each 4>k^x via heatbath accord- 
ing to 



-S^ff^bi'Ph,3i) 



m 



/j 0k,x-\-jj,^ 



-s{U'^^^,4>k,.,4>k,.+^) 



(11) 



and finally update (fik^x+p. with the standard heat- 
bath procedure. This set of moves is presently 
being implemented. 

To estimate the effects of this kind of update 
we simulated these moves by multiple heatbath 
updates on the fields Up^^x-, 4'k,x and (fik^x+p. and 



measured the probability distribution of the an- 
gle between the updated U' ^ and U^^^x ■ The re- 
sults are shown in Fig. 2. One sees that on one 
hand the width of the distribution is increased 
considerably, but that the n-dependence is still 
present. This means that at least the propor- 
tionality factor of n in rmt will be lowered con- 
siderably. This data already indicates that either 
the n-behaviour is changed or the proportionality 
constant is decreased by at least one to two orders 
of magnitude smaller than without modification. 
By implementing this move we expect to improve 
the efficiency of the algorithm considerably, even 
taking into account an increased computational 
cost. 

5. CONCLUSIONS 

While the systematic errors in the approxima- 
tions seem to be rather uncritical, as stated in 
[2], the Monte Carlo part of the algorithm needs 
to be improved further. Some ideas which will 
help considerably are being implemented, but es- 
pecially the local structure of the method makes 
the whole range of existing MC methods (except 
cluster algorithms) applicable. 
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